clear
clc
set(0,'DefaultFigureWindowStyle','docked')
addpath(genpath('library.model.20220926'))

%% about
% fits the following parameters:
% b, kappa, eta, z(1), 
% x_bar, sigma_x, chi
% to target the following moments:
% E[U], std[U], std[N*w]/std[Y], std[output growth vol], 
% E[phi|expansion], std[mu3_1y], std[1 year nom yld]

% in the data: E[s]=0.034 and std[s]=0.006
% Bond.Risk.Premia.Mitra.Xu\Data\RFS_RR_Empirics\all
% empirics\macro_moments_1964_2016.do

%% parameters

% grid sizes
parms.NN = 50;
parms.Nx = 50;
parms.Nz = 2;

% grid: N
N_lb = 0.85; 
N_ub = 0.999;
parms.grid_N = linspace(N_lb, N_ub, parms.NN)';

% grid: z
parms.Nz= 2;
parms.fract_bad = 1/6;
parms.rho_z = 0.9; % 0.9: recessions (booms) last 12 (60) months on average
parms.grid_z = [-0.03 0];

pi_LR(1) = parms.fract_bad;
pi_LR(2) = 1-pi_LR(1);
for s_prime = 1:2
    for s = 1:2
        parms.StateTransitionProbs(s,s_prime) = pi_LR(s_prime)*(1-parms.rho_z) + parms.rho_z*(s_prime==s);
    end
end

% labor market
parms.b = 0.81; % Unemployment benefit is b
parms.beta = 0.9979; %  Rep. household's time preference parameter 
parms.iota = 1.24; % Matching function
parms.phi_match = exp(-0.1/12); % constraint on matching (slow recovery)
% parms.phi_match = NaN; % constraint on matching (slow recovery), NaN=off
parms.kappa = 0.76; % vacancy costs
parms.eta = 0.158; % Bargaining power of worker

% parms.s = 0.034*ones(parms.Nz,1); % Hall AER (2017) following Shimer's calibration uses 0.0345

% parms.s = [0.0474;0.0313]; % E[s]=0.034, std[s]=0.006
parms.s = [0.037;0.033]; % see all_empirics.pdf

% preferences
parms.chi = 1/3.2; % CES aggregator
parms.gamma = 2; % risk aversion (intertemporal)

PHI_MAX = 0.99;
PHI_MIN = 0.1;
% PHI_BAR = 0.9;

phi_grid = linspace(PHI_MIN,PHI_MAX,parms.Nx)';

% parms.x_bar = -log(1/PHI_BAR - 1);
parms.x_bar = 2.5572;
% parms.rho_x = 0.5^(1/12);
parms.rho_x = 0.9^(1/12);
% parms.rho_x = 0.85^(1/12);
parms.grid_x = -log(1./phi_grid - 1);
parms.sigma_x = 0.13*ones(parms.NN,parms.Nz); % should have sigma>0

parms.OPTION_x_standardize_shock = true; % dx loads on z'-E[z']/(std(z'))
% parms.OPTION_x_standardize_shock = false; % dx loads on z'-E[z']

parms.sigma_x(:) = 0.1; % should have sigma>0

%% check parameters

CheckParams_Model_20220926(parms)

%% Fit model

NumTries = 1;
opts = optimset('display','iter-detailed','UseParallel',true,'MaxFunEvals',100);

U_mean_target = 0.0609;
U_std_target = 0.0078;
wagebillratio_target = 0.87;
vol_gy_target = 0.0081;
nom_yld_1y_mean_target = 5.28;
nom_yld_1y_vol_target = 3.32;
ave_mu3_vol_target = 0.0304*0.8^3;
E_phi_boom_target = 0.91;

kappa_guess = 0.15;
b_guess = 0.9;
eta_guess = 0.3;
z1_guess = -0.036;
beta_guess = 0.9979;
beta_lb = 0.995;
beta_ub = 0.999;
chi_guess = 1/3.2;
x_bar_guess = 2.5; 
log_sigma_x_guess = log(0.14);
    
% [parms, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
%     FitModel20220926_baseline(parms, kappa_guess, b_guess, eta_guess, z1_guess, chi_guess, x_bar_guess, log_sigma_x_guess, ...
%         U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_vol_target, E_phi_boom_target, ave_mu3_vol_target, NumTries, opts);

% [parms, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
%     FitModel20220926_baseline_ver2(parms, kappa_guess, b_guess, eta_guess, z1_guess, beta_guess, beta_lb, beta_ub, chi_guess, x_bar_guess, log_sigma_x_guess, ...
%         U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_mean_target, nom_yld_1y_vol_target, E_phi_boom_target, ave_mu3_vol_target, NumTries, opts);

% fix x_bar
[parms, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_baseline_ver3(parms, kappa_guess, b_guess, eta_guess, z1_guess, beta_guess, beta_lb, beta_ub, chi_guess, log_sigma_x_guess, ...
        U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_mean_target, nom_yld_1y_vol_target, ave_mu3_vol_target, NumTries, opts);

%% Verify solution
tic
Verify_Solution_Model_20220926(parms,gesol)
toc

%% compute model moments

[gesol,Phi_Nxz,prob_Nxz,Phi_NxYz,grid_Y,prob_NxYz,stats] = ...
    ComputeModelMoments_Model20220926(parms,gesol);

%% save

save([mfilename,'.mat'])